Transcriptome analysis using next generation sequencing reveals molecular signatures of diabetic retinopathy and efficacy of candidate drugs.

PURPOSE
To define gene expression changes associated with diabetic retinopathy in a mouse model using next generation sequencing, and to utilize transcriptome signatures to assess molecular pathways by which pharmacological agents inhibit diabetic retinopathy.


METHODS
We applied a high throughput RNA sequencing (RNA-seq) strategy using Illumina GAIIx to characterize the entire retinal transcriptome from nondiabetic and from streptozotocin-treated mice 32 weeks after induction of diabetes. Some of the diabetic mice were treated with inhibitors of receptor for advanced glycation endproducts (RAGE) and p38 mitogen activated protein (MAP) kinase, which have previously been shown to inhibit diabetic retinopathy in rodent models. The transcripts and alternatively spliced variants were determined in all experimental groups.


RESULTS
Next generation sequencing-based RNA-seq profiles provided comprehensive signatures of transcripts that are altered in early stages of diabetic retinopathy. These transcripts encoded proteins involved in distinct yet physiologically relevant disease-associated pathways such as inflammation, microvasculature formation, apoptosis, glucose metabolism, Wnt signaling, xenobiotic metabolism, and photoreceptor biology. Significant upregulation of crystallin transcripts was observed in diabetic animals, and the diabetes-induced upregulation of these transcripts was inhibited in diabetic animals treated with inhibitors of either RAGE or p38 MAP kinase. These two therapies also showed dissimilar regulation of some subsets of transcripts that included alternatively spliced versions of arrestin, neutral sphingomyelinase activation associated factor (Nsmaf), SH3-domain GRB2-like interacting protein 1 (Sgip1), and axin.


CONCLUSIONS
Diabetes alters many transcripts in the retina, and two therapies that inhibit the vascular pathology similarly inhibit a portion of these changes, pointing to possible molecular mechanisms for their beneficial effects. These therapies also changed the abundance of various alternatively spliced versions of signaling transcripts, suggesting a possible role of alternative splicing in disease etiology. Our studies clearly demonstrate RNA-seq as a comprehensive strategy for identifying disease-specific transcripts, and for determining comparative profiles of molecular changes mediated by candidate drugs.

Diabetes has emerged as a major worldwide public health concern, and the number of diabetics is estimated to exceed 400 million by the year 2030 [1]. A side effect of diabetes, namely diabetic retinopathy, is also a leading cause of blindness in working age adults (NIH MedlinePlus the Magazine). Several approaches, including good glycemic Correspondence to: Raj P. Kandpal control, use of blood pressure medications, and lipid control, have been demonstrated to inhibit diabetic retinopathy in clinical trials, but many patients are not able to maintain these regimens over the long-term. Thus, additional therapeutic approaches are continuously being sought. Several experimental therapies that include vitamin E, aspirin, aminoguanidine, or inhibitors of receptor for advanced glycation endproducts (RAGE) and p38 mitogen activated protein (MAP) kinase [2][3][4][5][6] have shown positive effects at inhibiting the development of diabetic retinopathy lesions in laboratory animals, but the underlying molecular mechanisms are not clear. Given their importance in cellular metabolism and regulatory processes, these therapeutic agents are expected to target distinct pathways either directly or indirectly. Therefore, identification of the targets of these drugs might assist in characterizing their molecular side effects.
Molecular changes accompanying the progression of disease can now be determined by several approaches. Gene expression microarray analysis has been widely used during the past decade for characterizing complete transcriptomes [7][8][9] and it has yielded global profiles of whole retina or retinal cell types in both wild type and disease models [10][11][12][13][14][15][16][17][18]. A comparison of the expressed complement of the genome between normal and diabetic retinas has indicated altered abundance of transcripts involved in several key pathways [19,20]. Although microarray strategies have been successful in describing disease-or phenotype-associated expression changes, hybridization-based profiling approaches suffer from technical variations that are difficult to control. For this reason, many expression changes cannot be validated by quantitative reverse transcription polymerase chain reaction (qRT-PCR) [21]. In addition, relevant causative expression changes, such as alternatively spliced variants of transcripts and the expression of novel transcripts in disease samples, may not be comprehensively captured because specific probe sets may not be included on the particular microarrays being used.
Next generation sequencing based on the RNA sequencing (RNA-seq) approach is now gaining prominence as a means of accurate qualitative and quantitative characterization of the expressed complement of a genome [22,23]. This method provides millions of sequences from expressed RNA molecules and can provide relatively unambiguous definition and abundance of transcripts in a given specimen. RNA-seq is therefore expected to reveal a better representation of the transcriptome, and this strategy is also more amenable for the analysis of alternatively spliced transcripts.
We have recently demonstrated the high accuracy and sensitivity of RNA-seq technology with microarray and qRT-PCR methods by profiling the neural retina specific leucine zipper deficient (Nrl −/− ) retina [24]. In this paper, we define the changes occurring in the entire retinal transcriptome of streptozotocin (STZ)-induced diabetic mice compared to wild type mice. We also attempt to identify relevant molecular signatures and cellular pathways, with the goal of examining the impact of inhibitors of RAGE and of p38 MAP kinase, two candidate drugs that inhibit retinopathy in experimental models. Our results indicate that these inhibitors induce common, as well as unique, molecular changes that may assist in determining the efficacy and/or side effects of candidate drugs.
METHODS Streptozotocin-induced diabetes in mice: All animal experiments conformed to the ARVO Resolution on the Use of Animals in Research and were approved by the Animal Care and Use Committee of Case Western University. Diabetes was induced in C57BL/6 mice (Jackson Laboratories, Bar Harbor, ME) with streptozotocin and a corresponding number of weight and age-matched animals were maintained as normal controls. Mice received five sequential daily intraperitoneal injections of a freshly prepared solution of streptozotocin (Sigma-Aldrich, St Louis, MO) in citrate buffer (pH 4.5) at 60 mg/kg of bodyweight. After hyperglycemia was verified at least three times during the second week after streptozotocin administration, diabetic mice were randomly assigned to remain as untreated diabetic controls or were administered RAGE-Ig fusion protein intraperitoneally at three different concentrations (10, 100, and 300 µg per mouse) three times per week. Insulin was given to diabetics as needed (0-0.3 units of neutral protamine Hagedorn (NPH) insulin subcutaneously, 0-3 times per week) to prevent weight loss without preventing hyperglycemia.
The RAGE fusion protein (provided by L. Brown, Galactica Pharmaceuticals, Inc., Villanova, PA) consists of a RAGE ligand binding element, a heavy chain immunoglobulin of G4 isotype (IgG4) constant domain with a linker connecting the ligand binding element with the constant domain. The RAGE-Fc binds to all of the known ligands of RAGE and acts as a competitive, negative regulator of RAGE signaling by competing with the membrane-bound receptor for binding of ligands. For the present studies, a murine version of hRAGE was employed in which a G2a isotype (IgG2a) constant domain (biologically equivalent to IgG4 in hRAGE) was used. The inhibitor of p38 MAPK (PHA666859) was provided by Pfizer Research Laboratories (Groton, CT). PHA666859 was mixed into powdered diet and replaced weekly, and food consumption was measured to calculate the amount of drug consumed.
We have previously shown that both of these therapies inhibit the development of early stages of diabetic retinopathy, as determined by the abundance of acellular capillaries and pericyte ghosts [3,25]. Glycated hemoglobin (GHb; an estimate of the average level of hyperglycemia over the previous 2-3 months) was measured by affinity chromatography (Glyc-Affin; Pierce, Rockford, IL) every 3 months in each animal, after an overnight fast . Bodyweight and average daily food consumption were measured weekly. Experimental variations were reduced by CO2 asphyxiation of all animals at 8 months of diabetes between 2 -4 PM Isolation of RNA: Retinal RNA was extracted with Trizol reagent [26,27]. Three independent retina samples were obtained from groups of animals that were normal, diabetic, treated with RAGE inhibitor, and treated with p38 MAPK inhibitor. The RNA was quantified by absorbance at 260 nm in a spectrophotometer and its integrity was determined using an Agilent Bioanalyzer (Santa Clara, CA) [28]. Illumina RNA-seq: Independent RNA samples were used for the preparation of cDNA libraries using an mRNA-seq Sample Preparation Kit (Illumina, San Diego, CA). Briefly, double stranded cDNA was prepared by random hexamer priming and reverse transcription of chemically sheared total RNA (2 µg). The double stranded cDNA was ligated with adaptor sequences required for bridge amplification using the Illumina flow-cell, and the ligated DNA was amplified before library assessment and quantitation using the BioAnalyzer DNA-1000 chip (Agilent). Cluster generation of 10 pM loads of a single mRNA-seq library per lane was performed on the Cluster Station using an Illumina Single Read Cluster Generation Kit v4 (Illumina). Sequence-by-synthesis single reads of 54 base length using the SBS Sequencing Kit v4 (Illumina) were generated on the Genome Analyzer IIx. This instrument runs real-time analysis SCS2.4 software (Illumina) for base calling. Real-time qPCR: Quantitative RT-PCR was performed using a Genetic Analyzer 7900HT (Life Technologies, Carlsbad, CA) with Taqman or SYBR green assays, as previously described [24]. Briefly, total RNA (1 µg) was reverse transcribed by Superscript II reagents (Life Technologies, Carlsbad, CA) primed with oligo(dT)20. The cDNA was diluted fivefold and one microliter was used as template in subsequent qRT-PCR reactions. The integrity of the PCR reaction was verified by melt curve analysis. The results of amplifications of three biologic replicates, performed in triplicate, were averaged to produce the qRT-PCR data. Changes in gene expression were determined by the ddCT method using the levels of HPRT (hypoxanthine guanine phosphoribosyl transferase) transcript for normalization. Bioinformatics analysis of RNA-seq data: The cDNA sequences captured on the Illumina platform were analyzed using the following workflows.
(1) Transcript isoform level analysis was performed by aligning the 54 base cDNA reads against the reference genome mm9 build using a Burrows-Wheeler transform based short read aligner (BWA) [29], and the aligned reads were visualized using Integrated Genome Viewer [30]. The read alignment files were imported into Partek Genomics Suite (Partek® Genomics SuiteTM. 6.3 ed. St. Louis, MO: Partek Inc.; 2008) and RPKM (reads per kilobase of exon model per million mapped reads) counts for each of the 28,157 transcripts defined in the UCSC refflat annotation file were calculated. A stringent filtering criterion with RPKM value of 1.0 [31] was used to obtain 16,577 expressed transcripts. The RPKM values of filtered transcripts were log-transformed using log2 (RPKM + offset) with an offset value of 1.0. This filtering protocol reduced the number of transcripts from 28,157 to 16,577 for further analyses. Fold changes in transcript expression and p-values were computed using ANOVA, and significantly altered transcripts were selected by applying fold change-cutoff of 1.5 and/or p-value cutoff of 0.05. Gene Ontology (GO) pathways that were significantly enriched between each of the sample pairs were identified using ExPlain tool from Biobase International, Beverly, MA.
(2) Alternatively, gene level analysis was performed by aligning the filtered reads to the mouse reference genome build mm9 using the ELAND (Efficient Local Alignment of Nucleotide Data) algorithm (Anthony J Cox, Solexa Ltd Saffron, Walden, UK). First, raw and RPKM normalized counts were calculated for gene models (as defined in UCSC RefGene). Subsequently, we applied SAM (significance analysis of microarrays) to gene level RPKM to identify RNAs with a fold change greater than 1.5 and false discovery rate (FDR) of less than 25%. Briefly, the boundaries of exons were obtained from the RefGene database (Genome), and the numbers of mapped reads on each exon were calculated. The mapped reads as RPKM indicated the expression level of exons, and this value was used to calculate the Log base 2 of the average (or maximum) expression level of exons for each gene. For this alternative analysis, a less stringent cutoff of 0.5 RPKM was used to select expressed genes for further analysis. Alternatively spliced exons: Alternatively spliced transcripts were detected using the AltAnalyze tool [32], which employs de novo splice junctions to predict alternative splicing events based on the reads spanning the splice junctions. First, TopHat [33] was run independently on sequence reads obtained for each of the 12 RNA samples to identify all novel splice junctions. Second, we combined all predicted splice junctions into a master file and used this as an input option, together with the Ensembl mouse gene transfer format (GTF) annotation file, for a second iteration of TopHat runs on each of the 12 samples. Finally, these 12 junction files were used as input in AltAnalyze to detect alternative splicing using the analysis of splicing by isoform reciprocity (ASPIRE) algorithm.
Briefly, we calculated the number of sequences mapping to a specific exon junction and the numbers of sequences mapping to junctions as well as exons of that specific gene. These two numbers were then used to calculate the fraction of the junction containing sequences relative to the mean of abundance of sequences of all junctions and exons in the gene. The procedure was applied to each reciprocal isoform pair (i.e., an isoform that includes a specific exon and an isoform in which that exon is excluded). The inclusion and exclusion ratios were then used to compute the value for ΔI. The details of these protocols have been described in AltAnalyze workflow [32], as adapted from ASPIRE algorithm [34]. The ratios described above were used to calculate p values for the sample ASPIRE scores relative to controls, and to determine false-discovery rate p values [35].

Glycemia:
Diabetic mice from all experimental groups had levels of GHb and blood glucose that were significantly greater (p<0.05) than levels found in appropriately agematched nondiabetic controls. Average GHb for the nondiabetic control (N), diabetic control (D), diabetic plus PHA666859 (D + PHA666859), and diabetic plus murine RAGE-Fc fusion protein (D + mRAGE-Fc) groups over the entire duration of study are listed in Table 1. Although diabetic mice did not gain weight at the normal rate, none of the animals lost bodyweight, and all of the animals appeared clinically healthy. These results indicate that the therapies involving inhibitors of RAGE and p38 MAP kinase do not alter the diabetic status of the animals. Raw read counts and transcripts in sequenced samples: Total raw reads for various samples were in the range of 27.8-35.5 million, which represented ~1.6 gigabases of sequences for each sample. Efficient Local Alignment of Nucleotide Data (ELAND) and Burrows-Wheeler Alignment (BWA) tools allowed the alignment of 13.8-20.9×10 6 reads from various samples to mouse genome sequence database (build mm9). These aligned reads were subsequently analyzed by significance analysis of microarrays (SAM) or ANOVA. The sequence alignment files were imported into PARTEK, and the normalized values for numbers of RPKM for 28,157 transcripts defined in the UCSC refflat file were obtained. Among these sequences obtained from nondiabetic and diabetic animals, 16,231 transcripts had RPKM values of greater than one, and included rhodopsin (Rho), guanine nucleotide binding protein (Gnat1), retinol binding protein (Rbp1), phosphodiesterase 6A (Pde6a), and guanine nucleotide binding protein (Gngt1) as the most abundant sequences specific to photoreceptor cells. Similarly, the list also contained other less abundant transcripts such as pairedlike homeodomain 3 (Pitx3), PDZ and LIM domain 1(Pdlim1), and meteorin (Metrnl). The presence of these highly abundant as well as rare transcripts in the data set demonstrates that the library preparation and data analysis methods did not introduce any undesirable bias. Comparison of transcript differences revealed by two analyses: We analyzed the sequence data using ELAND-SAM and BWA-ANOVA methods, and compared the results to confirm method-specific changes in transcript abundance in the various groups of animals. Differential expression in diabetic and normal animals revealed the top 100 transcripts shortlisted by the two methods, as shown in Table 2 (Appendix 1). Of these 100 transcripts, 81 are present in both data sets ( Figure 1) including a variety of crystallin genes. However, UDP-glucuronosyl transferase transcripts were detected only by the ELAND-SAM method. A significant difference was also noted in the fold-changes indicated by the two methods. These differences can be attributed to differences in the alignment algorithms, filtering criteria, and calculation of differential gene expression between the two methods, and demonstrate that altered abundance of transcripts can be detected by either one of these two methods.
Validation of significantly altered transcripts: Some representative transcripts were chosen for validation based on the magnitude of the changes in their abundance and their involvement in various pathways. We performed Taqman real-time PCR to validate 11 transcripts, including Crygacrystallin gamma A, Crygb-crystallin gamma B, Dctdopachrome tautomerase, Egr-early growth response, FABPfatty acid binding protein, Lgsn-lengsin, Lim2-lens intrinsic membrane protein, Nr2c2ap-nuclear receptor 2C2-associated protein, retinal pigment epithelium 65 Kd protein (Rpe65), Sfrp1-secreted frizzled-related protein, and wingless type MMTV integration site, member 7b (Wnt7b). The level of each of these transcripts corresponded to the pattern revealed by RNA-seq, and indicated significant changes in the levels of these transcripts between diabetic and normal mice ( Figure  2). However, the magnitude of change observed by RT-PCR was not identical to the results of RNA-seq. The noncorrespondence of magnitude reflects the differences in sensitivities of the RNA-seq and Taqman assays, which are two inherently different strategies.
A few additional transcripts (Cdh3-cadherin, Fgf2fibroblast growth factor, Lrat-lecithin retinol acyltransferase, Mnd1-meiotic nuclear divisions 1 homolog, Rgr-retinal Gprotein coupled receptor, and Sema3c-semaphorin) were chosen for validation with SYBR green based real-time qPCR. All of these transcripts were significantly altered in diabetic mice ( Figure 3). As stated for Taqman qPCR validation, the correspondence between RNA-seq and SYBR green qPCR was comparable but not identical.
Transcripts, pathways and gene ontology (GO) categories are altered in retinas isolated from diabetic animals: The application of SAM and a false discovery rate of 25% for the data set revealed greater than 1.5 fold alterations in the levels of 173 transcripts. Alternative analysis with BWA-ANOVA, however, indicated alterations of greater than 1.5 fold in 298 *The values are reported as mean ± standard deviation, and were calculated for 10, 7, 7, and 8 animals in the four groups, respectively. non-diabetic controls. These transcript changes were detected by both analysis methods.
Alterations in Wnt signaling transcripts: Altered levels were observed for 14 transcripts that are either directly involved in Wnt signaling pathway or are targets of Wnt signaling (Table  4). While Sfrp1 and Wnt7b levels in diabetic animals increased by 3.9 and 2.9 fold, Wnt7a and Mfrp were upregulated by more than 2.0 fold. The increase in the levels of Wnt5b and axin, and decrease in Wisp 1 and Wnt2b transcripts were marginal. The levels of a few target genes of Wnt signaling such as FRAT2, BMP4, Cldn1, and Cldn2 were also altered.
Transcripts involved in microvasculature and the inflammatory pathway: Endothelin and VEGF were among the vascular transcripts detected in the retina. The levels of Edn2 decreased by 1.4 fold, whereas an increase of 1.4 fold was observed in the levels of Edn3 (Table 5). The levels of other vasculature related genes were marginally altered in diabetic animals. Alterations were also observed in the abundance of inflammation associated transcripts, namely Ltbp1, Bmp4, Hspb1, CD44, C1qtnf5, Ifitm1, and Islr, which varied between 1.5 and 2.2 fold in diabetic animals ( Table 6).

Transcripts involved in apoptosis:
We analyzed a variety of pro-apoptotic and anti-apoptotic transcripts as well as other transcripts that influence the apoptosis pathways. A significant increase was observed in the levels of Dapl1 and Bid1, while the levels of Fads3, Fas, and Bag5 showed marginal alterations. An appreciable decrease was also noted in the levels of Trafip3, Tnfsf13, Tnfrsf18, and Casp7 (Table   7). It warrants mention that these transcripts were not shortlisted when the FDR and p-value cut-off were applied to the data analysis. Given the importance of apoptosis in diabetic retinopathy [36,37], these genes were specifically analyzed from the master list of sequenced tags and were found to have altered levels in the diabetic animals.
Transcripts involved in neuronal functions and photoreceptors: The list of differentially expressed genes obtained by applying stringent threshold parameters did not contain any transcripts involved in neuronal development or signal transduction. However, a review of the transcript master list, without applying rigorous statistical significance, revealed changes in transcripts encoded by polycomb genes (Pcgf), acetyl choline receptor genes (Chrn), muscarinic cholinergic receptor genes (Chrm), potassium channels (Kcnq), glutamate receptor genes (Grin), solute carrier/ transporter genes (Slc), and the fasciculation and elongation gene (Fez) ( Table 8). Similar review revealed decreased levels of transcripts involved in the development and function of photoreceptors. These transcripts included cyclic nucleotide gated channel (Cngb3), arrestin (Arr), guanine nucleotide binding protein (Gnb3), and phosphodiesterase (Pde6h). Marginal decreases were also noticed in photoreceptorspecific transducin (Gnat2) and Crxos1. However, the levels of Nrl, Crx, and Nr2e3 remained relatively unchanged (Table  9).
Altered expression of genes encoding glucuronosyl transferases in diabetic retinopathy: A significant increase Figure 3. Cyanine (SYBR) green realtime quantitative polymerase chain reaction (qPCR) validation was performed for selected transcripts. The transcripts were amplified as described in Methods section. Three technical replicates were performed for three biologic replicate samples of RNA isolated from nondiabetic animals, diabetic animals and diabetic animals treated with either inhibitor of receptor for advanced glycation endproducts (RAGE) indicated as RI, or an inhibitor of p38 mitogen activated protein kinase (MAPK) designated as p38. The bars (left to right) represent samples corresponding to diabetic animals, diabetic animals treated with RAGE inhibitor, and diabetic animals treated with p38 MAPK inhibitor, respectively. The error bars represent standard deviation from the mean value (n=9). The levels of transcript in diabetic and treated animals were determined relative to the levels present in nondiabetic controls.       was observed in transcripts for members of the UDPglucuronosyl transferase gene family during diabetic retinopathy. The protein products of these genes are primarily responsible for the metabolism of xenobiotics. In particular, UGT1A5, UGT1A10, UGT1A1, UGT1A9, UGT1A2, UGT1A7C, UGT1A6A, and UGT1A6B transcripts were elevated by greater than 3.0 fold in diabetic animals (Table  10). These transcripts were detected in analyses performed by the ELAND-SAM method at the gene level with a FDR cutoff of 25%, but were absent in the list of transcripts generated by transcript isoform level BWA-ANOVA analysis at a p value cutoff of less than 0.05.

Rhodopsin-retinol pathway altered in diabetic mice:
Rhodopsin-retinol pathway transcripts showed increased levels for lecithin retinol acetyl transferasease, transthyretin, retinol binding protein, and Stra6 (stimulated by retinoic acid). However, decreases occurred in the levels of transcripts corresponding to retinoic acid metabolism enzyme Cyp26a1, retinol dehydrogenase, and rhodopsin phosphorylating GPCR kinase (Grk4).

Pathways related to altered metabolism in diabetic mice:
Transcripts directly/indirectly related to glycemic pathway showed marginal increases in the case of glycerol-3phosphate dehydrogenase and insulin secretion modulating Trpm4 protein, whereas transcripts levels decreased slightly for a glycogen interacting protein Trim7 and for insulin-like growth factor binding proteins Igfbp3 and Igfbp5. Alternative splicing of genes in diabetic retinopathy: The RNA-seq method provides the identity of transcripts that have yielded an unambiguous sequence. Some of the sequence reads, which span exon boundaries, can be analyzed to predict alternatively spliced transcripts. We identified 55 genes that were alternatively spliced in the retina of diabetic animals ( Table 11). The alternative splicing events in these genes were detected by selecting transcripts with a RPKM cutoff of 3.0 and using ASPIRE's ΔI cutoff of 0.2, which indicates a 20% change in exon inclusion. This threshold of ΔI has been shown to yield accurate representations of alternatively spliced versions of transcripts [34].
It warrants mention that alternative splicing of transcripts is also prevalent in the retinas of normal animals. The 55 alternatively spliced transcripts therefore represent either novel spliced versions or spliced isoforms with altered abundance in diabetic animals, which suggests an association of alternative splicing with diabetic retinopathy. Noteworthy among these transcripts were those involved in neuronal cell survival, neural differentiation, apoptosis, endocytosis, intracellular transport, glucose metabolism and insulin impairment, G-protein coupled receptor signaling, and mRNA processing. Specifically, alternative splicing of arrestin (Arr), interphotoreceptor matrix proteoglycan (Impg2), and transient receptor membrane potential cation channel provide tentative evidence for the effects of diabetic retinopathy on the activity and integrity of photoreceptor cells in the retina.

Responsiveness to treatments and molecular side effects of therapeutic intervention with candidate drugs:
We have characterized the molecular responsiveness of animals to treatments as positive, negative, or as side effects, based on a comparison with wild type mice. We define positive effects of drugs as a reversal of transcript levels toward normalcy, and negative effects as changes in transcript levels of treated animals that are greater than those seen in the diabetic animals when compared to wild type animals. If a transcript remains unchanged in retinas from diabetic animals and the therapy leads to a significant change in its abundance, this is classified as a side effect.
Treatments with inhibitors of RAGE and p38 MAPK inhibited diabetes-induced changes (positive effects) and also exacerbated other changes (negative effects; Table 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, and Table 10). Both inhibitors blocked the diabetes-induced changes in levels of all crystallin transcripts and the levels of UDPglucuronosyl transferases. Both inhibitors worked well in preventing the changes in the levels of Wnt pathway transcripts or apoptosis transcripts, but the RAGE-inhibitor appeared to work positively on a greater number of transcripts, as was the case for transcripts of inflammatory pathways. Surprisingly, photoreceptor transcripts were adversely affected in animals treated with either of these drugs. Of the 837 transcripts that were altered by greater than 1.5 fold in diabetic mice or in any one of the treatments, 379 transcripts were further exacerbated by RAGE inhibitor. The treatment with p38 MAPK inhibitor, on the other hand, caused further adverse exacerbations in 535 transcripts. Whether these exacerbations of transcript levels by drugs were beneficial or detrimental is not known at present. Similar analyses were also performed with the data related to alternative splicing ( Figure 4). Comparison of diabetic versus nondiabetic, RAGE inhibitor-treated diabetic versus nondiabetic and p38 MAPK inhibitor-treated diabetic versus nondiabetic animals led to the identification of alternatively spliced transcripts specific to each group of animals (Table 11, Table 12, and Appendix 2). Some notable alternatively spliced transcripts in diabetic animals included doublecortin-like kinase Dclk3, death effector domain Dedd2, dynein Dync1li2, interphotoreceptor matrix protein Impg2, ankyrin 2, arrestin 3, coiled-coil domain containing Ccdc64, G-protein coupled receptor Gprasp2, Nsmaf, dynactin 6, and whirlin (Table 11). Among the transcripts responsive to RAGE inhibitor were photoreceptor specific gene arrestin, Nsmaf transcript involved in apoptosis, and endocytosis transcript Sgip1. The exacerbated levels of 76 alternatively spliced transcripts in p38 MAPK inhibitortreated animals included the Wnt pathway transcript Axin1 and its target gene Cldn1, opsin transcript Opnm1w, Mapk8ip3, peroxisomal transcript Pex16, cilia transcript Cep164, and ocular regulator Pax6 (Table 12 and Appendix 2). A comparison of alternatively spliced transcripts revealed 12 transcripts common among nondiabetic animals, diabetic animals, diabetic animals treated with RAGE inhibitor, and diabetic animals treated with p38 MAPK inhibitor (Figure 4). Among others, these genes included G-protein coupled receptor associated sorting protein, kinesin family member 1A, nuclear receptor co-repressor, phospholipase A2, and zinc finger protein 444. Although the summary presented in Figure  4 does not allow us to decipher the mechanisms of these therapies, the results suggest that the changes in splicing patterns in response to therapy may reflect the positions of p38 MAP kinase and RAGE in the cell signaling pathways relative to each other.

DISCUSSION
This is the first report to describe the application of RNA-seq for comprehensive sequencing of transcripts in a diabetic retinopathy model and for evaluating the efficacy of candidate drugs. The RNA-seq methodology has allowed accurate and quantitative identification of molecular signatures. The quantitative RT-PCR data (by TaqMan or SYBR Green) validated the expression changes for all 15 transcripts that were examined. Interestingly, the rank order of differentially    expressed transcripts (as a consequence of diabetes or treatment with drugs) differed between the two analysis methods used. The ELAND-SAM method used gene model annotations (RefGene) from UCSC, while the BWA-ANOVA method used the transcript model annotations (refflat) from UCSC. However, the transcripts revealed by both methods of analyses were verifiable by real-time PCR and thus the validity of both methods was supported. The validation of the RNA-seq platform for comprehensive transcriptome analysis has recently been described by Brooks et al. [24]. Given the growing use of next generation sequencing platforms, the development of additional analytical tools in the future should allow even better description of RNA-seq data.
As predicted, we observed changes in a variety of retinal transcripts due to diabetes. Some of these changes have been found to influence glucose metabolism or to alter the levels of hormone and hormone receptors in other tissues; these include glycerol-3-phosphate dehydrogenase, galanin receptor (Galr2) involved in G-protein coupled receptor signaling, transient receptor (Trpm4) involved in insulin secretion, insulin like growth factor binding proteins (Igfbp3 and Igfbp5), the kinesin family of proteins, and calpain. Since vascular cells occupy such a small portion of the retina, the majority of observed changes likely are due to changes in the retinal neuroglia. Consistent with this idea, the list of transcripts that were significantly altered in diabetes did not contain transcripts corresponding to retinal vasculature, but targeted analyses of specific transcripts (such as endothelin-2, endothelin-3 and VEGFB) nevertheless revealed diabetesinduced alterations of 1.1-1.4 fold in levels of these transcripts.
Inflammation has been previously demonstrated to play an important role in the pathogenesis of diabetic retinopathy in animals [38,39]. Previous microarray profiling in a rat diabetic retinopathy model indicated changes in the levels of a variety of transcripts involved in inflammatory pathways [19]. Our RNA-seq data showed many inflammationassociated genes with altered expression; these included C1q TNF-related gene 5 (C1qtnf5), latent TGF-β binding protein (Ltbp1), phospholipase (Pla2g2f), interleukin 17 receptor 7c (Il17rc), and C4b. The changes in relatively few transcripts of the inflammatory pathway by eight months of diabetes seen in our study, when compared to changes reported at shorter durations, might indicate that inflammation is most severe early during the etiology of diabetic retinopathy, and that its severity might diminish as the disease progresses.
A variety of transcripts were altered in late stage diabetes; our RNA-seq data identified the Wnt signaling pathway, crystallins, and various transcription factors. Significant alterations were observed in the abundance of transcripts corresponding to Wnt7b, secreted frizzled receptor related protein-1 (Sfrp1), dapper antagonist of β-catenin (Dact1), Wnt inducible secreted protein-1 (Wisp1), and membrane frizzled related protein (Mfrp). Modulation of the Wnt pathway has been reported as a likely target for intervention for diabetic retinopathy [40], and Wisp1 is known to upregulate the expression of the anti-apoptotic Bcl-X(L) gene [41]. In addition to the involvement of Wnt pathway proteins in other inflammatory diseases [42], alterations in Wnt signaling proteins have been shown in a rat diabetic retinopathy model [43].
We and others have previously observed a significant upregulation of different crystallin transcripts in diabetic retina [44]. The levels of fourteen crystallins increased significantly in diabetic mice. These proteins were initially characterized as structural components of the eye lens [45,46], but expression of crystallins has been shown in other tissues including the retina [47]. Crystallins that belong to the small heat shock protein family have a non-structural role [48], and retina-specific elevated levels of Cryaa,Cryba1,Crygs,Crybb2,Crygb, and Crygc have been reported in experimental uveitis [49]. Cryaa was shown to localize to photoreceptor cells in experimental uveitis, which likely represents a function in these cells to protect against apoptosis [49]. Our results clearly indicate a dramatic upregulation of various crystallins in diabetic retinopathy and suggest this to be a stress response, which may be similar to the reported effects of increased levels of Hspb1 transcript in retina [50]. Post-translational modifications have been reported to alter chaperone activity of crystallins [51][52][53], while glycation associated with aggregation is likely to compromise the protective effects of crystallins in diabetic animals. Transcriptional upregulation of crystallins may therefore possibly represent a stress response, but the translated proteins may be rendered dysfunctional by diabetes-induced and deleterious post-translational modifications and other crosslinking events.
The deleterious changes to the retina that occur in diabetes are likely to be associated with alterations in the levels of transcripts associated with apoptosis, survival, and metabolic injury. Consistent with this hypothesis, we have observed changes in a variety of transcripts associated with apoptosis pathways. The metabolic injury to the neural retina is suggested by the altered levels of transcripts that modulate oxidative stress in cells. Changes in transcripts such as cytochrome c oxidase, glutathione-S-transferase omega-1 and omega-2, microsomal glutathione S-transferase, and glutathione synthase indicate that cellular injury to the retina occurs due to oxidative stress. While cytochrome oxidase is a sink for molecular oxygen channeled through electron transport chain in mitochondria, glutathione synthase and glutathione-S-transferases are important components of redox reactions. The comparative analysis of gene level expression data also revealed differential expression of UDP-glucuronosyl transferases in diabetic animals. The UDP-glucuronosyl transferases metabolize a variety of xenobiotics, but they are also important for reactions occurring in cellular metabolism, such as those involving steroid hormones. We speculate that these enzymes may possess some other signaling activities, in line with the reported observation that UDP-glucuronosyl transferase 2B17 depletion led to downregulation of the Mcl-1 gene [54]. If these types of signaling activities associated with UDP-glucuronosyl transferases are involved in the transcription of diabetic retinopathy-related genes, then targeted inhibition of specific UDP-glucuronosyl transferases might have therapeutic implications for diabetic retinopathy.
Numerous groups have been exploring ways to inhibit diabetic retinopathy by blocking specific metabolic pathways. Although none of these targeted approaches have been successful in the clinic to date, several candidate drugs are being evaluated for efficacy. The alterations in the levels of superoxide, nitric oxide, cyclooxygenase-2, VEGF, and NFκβ in the early stages of diabetic retinopathy are consistent with the inflammation observed in retinal tissue. Our drug studies suggest that diabetes-induced alterations observed in inflammatory pathway transcripts are mediated in particular via p38 MAPK [40,55,56]. Another metabolic consequence of chronic hyperglycemia, in addition to inflammation, is the accumulation of advanced glycation end products (AGEs) [57]. These AGEs and other ligands are recognized by RAGE, and activated RAGE is associated with a variety of downstream consequences such as upregulation of VEGF and NF-κβ, and increased leucocyte adhesion in retinal microvascular endothelial cells, all of which are consistent with inflammation. Importantly, inhibitors of RAGE and p38 MAPK have been shown to inhibit the development of early stages of diabetic retinopathy in animals [3,25].
Since the phenotypic benefits of the two therapies are similar, assessing their efficacies against the molecular changes that they cause might provide important information. We have tested these effects to determine whether the drug prevents diabetes-induced changes in transcripts or further exacerbates diabetes-induced changes in transcript levels. Our data suggest that the RAGE inhibitor reversed the levels of a large number of transcripts whose levels were significantly altered in diabetic mice. Although the p38 MAPK inhibitor seemed comparably effective at inhibiting the development of diabetic retinopathy, it did not inhibit as many diabetesinduced changes in levels of transcripts, in contrast to what was seen with the RAGE inhibitor. One possible interpretation of this difference is that the transcripts normalized by inhibition of RAGE, but not p38 MAPK, are not critical for the development of the lesions that characterize early diabetic retinopathy. Interestingly, the p38 MAPK inhibitor caused further exacerbations in many more diabetes-induced changes in transcript levels than did the RAGE inhibitor. The intracellular protein p38 MAPK is downstream of RAGE [58] in one of the three pathways, such as PI3K, NF-kB, and MAPK that are initiated by AGE. Thus, one might expect a greater number of altered transcripts upon activation of RAGE, and reversal of transcript changes in animals treated with the inhibitor of RAGE.
An important and useful outcome of the RNA-seq approach was the identification of alternatively spliced transcripts during the development of diabetic retinopathy. Disease specific alternatively spliced transcript variants are shown to exist in various conditions [58,59], and our experimental studies show that alternatively spliced transcripts arise as diabetic retinopathy develops. Characterization of alternatively spliced transcripts that are specifically present or absent in a specific condition provides another opportunity for therapeutic intervention of diabetic retinopathy. Tocotrienol and EGCG are known to target splice sites and to modulate the abundance of alternatively spliced transcripts [60,61]. Exploration of the ability of these agents to inhibit diabetic retinopathy will be an interesting avenue of further research.
Reversal of the levels of alternatively spliced exons upon treatments also provides an additional validation of drug efficacy. Our novel finding of alternative splicing of transcripts in diabetes raises the possibility that these changes contribute to the pathogenesis of the retinopathy. From this perspective, the identification of photoreceptor-specific changes in arrestin, Impg2, and Trpm1 raises the further possibility of a contribution by photoreceptors to diabetic retinopathy. These possibilities remain to be tested experimentally.
Just as therapies were shown above to promote further exacerbation of some diabetes-induced changes in transcript levels, therapy-specific exacerbation of alternative splicing also occurs. These side effects were more prominent following treatment with the p38 MAPK inhibitor than with the RAGE inhibitor, but the significance of these differences is not known at present. It warrants mention that alternatively spliced transcripts have been identified in a variety of diseases [62], and that a specific spliced form of a gene such as COX-1 is one target of acetaminophen [59].
In conclusion, we have applied RNA-seq strategy and two analysis methods to determine the molecular signatures that involve representative transcripts of retinal vasculature, inflammatory pathway, Wnt signaling, apoptosis, crystallins, and UDP-glucuronosyl transferases. The importance of several of these abnormalities in the development of diabetic retinopathy has been identified previously, illustrating the power of the RNA-seq technique to identify important changes in disease. Of importance, however, is the identification of disease and treatment-specific alternative splicing by this strategy. The remarkable longevity of changes observed in transcript levels during the 8 months of diabetes in our study suggests that these changes may be necessary to promote and maintain the complications of diabetes. The inhibitors of RAGE and p38 MAPK reversed some of the diabetes-induced levels of various transcripts to normal levels, and exacerbated the changes in others. The significance and mechanisms underlying these alterations remain to be determined, but our results suggest RNA-seq as a desirable strategy for investigating the pathogenesis of diabetic retinopathy or other diseases, to distinguish between candidate drugs and to determine their mechanisms of action.